Entrega II

Author

C. Tangana - DNI: 00000000-X

Instrucciones (leer antes de empezar)

  • Modifica de la cabecera del documento .qmd tus datos personales (nombre y DNI). IMPORTANTE: no modifiques nada más de la cabecera.

  • Asegúrate, ANTES de seguir editando el documento, que el archivo .qmd se renderiza correctamente y se genera el .html correspondiente en tu carpeta local de tu ordenador (pulsa el botón Render o guarda el documento con Render on Save activado).

  • Los chunks (cajas de código) creados están o vacíos o incompletos, de ahí que la mayoría tengan la opción #| eval: false. Una vez que edites lo que consideres, debes ir cambiando cada chunck a #| eval: true (o quitarlo directamente) para que se ejecuten.

  • Recuerda que puedes ejecutar chunk a chunk con el botón play o ejecutar todos los chunk hasta uno dado (con el botón a la izquierda del anterior).

  • IMPORTANTE: solo se podrá subir un archivo .html al campus, no se evaluarán entregas sin el .html correctamente generado. Recuerda: el código es una herramienta, no el fin en esta asignatura. Se evaluará especialmente las interpretaciones y conclusiones que detalles. Un ejercicio con el código perfecto pero sin ningún tipo de razonamiento, interpretación o conclusión no superará el 4 sobre 10 de nota.

Paquetes necesarios

Necesitaremos los siguientes paquetes (haz play en el chunk para que se carguen):

rm(list = ls()) # Borramos variables de environment

# descomentar si es la primera vez (y requieren instalación)
# install.packages("tidyverse")
# install.packages("performance")
# install.packages("olsrr")
# install.packages("corrr")
# install.packages("corrplot")
# install.packages("skimr")
library(tidyverse)
library(rsample)
library(performance)
library(olsrr)
library(corrr)
library(corrplot)
library(skimr)

Caso práctico: análisis de datos de cáncer

El archivo de datos a usar lo cargaremos desde el csv cancer.csv.

# no cambies código
datos <- read_csv(file = "./cancer.csv")
datos
# A tibble: 3,047 × 13
   deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
       <dbl>     <dbl>      <dbl>          <dbl>       <dbl>     <dbl>
 1      165.     61898     260131           11.2       500.       39.3
 2      161.     48127      43269           18.6        23.1      33  
 3      175.     49348      21026           14.6        47.6      45  
 4      195.     44243      75882           17.1       343.       42.8
 5      144.     49955      10321           12.5         0        48.3
 6      176      52313      61023           15.6       180.       45.4
 7      176.     37782      41516           23.2         0        42.6
 8      184.     40189      20848           17.8         0        51.7
 9      190.     42579      13088           22.3         0        49.3
10      178.     60397     843954           13.1       428.       35.8
# ℹ 3,037 more rows
# ℹ 7 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
#   PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
#   BirthRate <dbl>, region <chr>

Los datos representan diferentes características socioeconómicas de distintas regiones, extraídas de the American Community Survey (census.gov), clinicaltrials.gov, y cancer.gov.

¿Nuestro objetivo? Ser capaces de predecir de manera lineal y MULTIVARIANTE la variable deathRate, nuestra variable objetivo, que representa la mortalidad media de cancer, por cada 100 000 habitantes. El resto de variables son:

  • medianIncome: mediana de los ingresos de la región.
  • popEst2015: población de la región
  • povertyPercent: porcentaje de población en situación de pobreza.
  • studyPerCap: ensayos clínicos relacionados por el cáncer realizados por cada 100 000 habitantes.
  • MedianAge: edad mediana.
  • region: nombre de la región.
  • AvgHouseholdSize: tamaño medio de los hogares.
  • PercentMarried: porcentaje de habitantes casados.
  • PctHS25_Over: porcentaje de residentes por encima de los 25 años con (máximo) título de bachillerato.
  • PctUnemployed16_Over: porcentaje de residentes de más de 16 años en paro.
  • PctPrivateCoverage: porcentaje de residentes con seguro de salud privado.
  • BirthRate: tasa de natalidad.

IMPORTANTE: siempre que puedas, relaciona los valores numéricos de la salida de R con su fórmula.

Ejercicio 1 (1.25 puntos)

Chequea que no hay problemas de rango/codificación: en caso de que alguna variable tenga dichos problemas, debes sustituir los valores fueran del rango razonable por su media o mediana (sin esos datos fuera del rango), justificando por qué usas la media o por qué usas la mediana, ejecutando el código que consideres (si te atascas, sigue haciendo con el dataset original)

Como se vio en la entrega anterior, por los rangos aportados, el único problema de rango/codificación lo tenemos en MedianAge, con valores por encima de un rango razonable.

datos |> skim()
Data summary
Name datos
Number of rows 3047
Number of columns 13
_______________________
Column type frequency:
character 1
numeric 12
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
region 0 1 4 9 0 4 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
deathRate 0 1 178.66 27.75 59.70 161.20 178.10 195.20 362.80 ▁▇▆▁▁
medIncome 0 1 47063.28 12040.09 22640.00 38882.50 45207.00 52492.00 125635.00 ▇▇▁▁▁
popEst2015 0 1 102637.37 329059.22 827.00 11684.00 26643.00 68671.00 10170292.00 ▇▁▁▁▁
povertyPercent 0 1 16.88 6.41 3.20 12.15 15.90 20.40 47.40 ▃▇▃▁▁
studyPerCap 0 1 155.40 529.63 0.00 0.00 0.00 83.65 9762.31 ▇▁▁▁▁
MedianAge 0 1 45.27 45.30 22.30 37.70 41.00 44.00 624.00 ▇▁▁▁▁
AvgHouseholdSize 0 1 2.48 0.43 0.02 2.37 2.50 2.63 3.97 ▁▁▃▇▁
PercentMarried 0 1 51.77 6.90 23.10 47.75 52.40 56.40 72.50 ▁▂▇▇▁
PctHS25_Over 0 1 34.80 7.03 7.50 30.40 35.30 39.65 54.80 ▁▂▇▇▁
PctUnemployed16_Over 0 1 7.85 3.45 0.40 5.50 7.60 9.70 29.40 ▅▇▁▁▁
PctPrivateCoverage 0 1 64.35 10.65 22.30 57.20 65.10 72.10 92.30 ▁▂▇▇▂
BirthRate 0 1 5.64 1.99 0.00 4.52 5.38 6.49 21.33 ▂▇▁▁▁
ggplot(datos) + 
  geom_histogram(aes(x = MedianAge)) +
  theme_minimal()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

En el histograma vemos que hay un corte muy claro, por lo que estableceremos (por ejemplo) un máximo de 120 años: por encima de dicha cantidad, se considerará una observación mal codificada. En lugar de eliminarlas lo que vamos a hacer es sustituirla por una medida de centralidad (imputación), que podemos hacer en base a dos argumentos:

  • Si usamos todos los datos, el histograma anterior, tenemos una distribución muy asimétrica, por lo que deberíamos usar la mediana (que se ve menos afectada por esos valores mal codificados, ya que es más robusta)
ggplot(datos) + 
  geom_density(aes(x = MedianAge)) +
  theme_minimal()

ggplot(datos) + 
  geom_boxplot(aes(x = MedianAge)) +
  theme_minimal()

  • En esta solución usaremos la media PERO sin esos datos (ya que en realidad…deberían ser NA ya que están mal recopilados).
ggplot(datos |> filter(MedianAge < 120)) + 
  geom_histogram(aes(x = MedianAge)) +
  theme_minimal()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(datos |> filter(MedianAge < 120)) + 
  geom_density(aes(x = MedianAge)) +
  theme_minimal()

ggplot(datos |> filter(MedianAge < 120)) + 
  geom_boxplot(aes(x = MedianAge)) +
  theme_minimal()

Vemos en los gráficos que es bastante simétrica y, de hecho, su distribución se asemeja mucho a una normal, por lo que la media (sin esos datos) si es representativa de la población, y lo que haremos será sustituir los valores por encima del límite permitido por dicha media

media <-
  datos |>
  filter(MedianAge < 120) |>
  summarise(mean = mean(MedianAge)) |>
  pull(mean)

# Escribe el código que consideres
datos_preproc <-
  datos |>
  mutate(MedianAge = if_else(MedianAge < 120, MedianAge,
                             media))

Ejercicio 2 (0.75 puntos)

En caso de encontrar alguna variable cualitativa, procesa los datos como consideras para incluir su información de manera que pueda ser usada en la futura regresión. Recuerda: si tenemos una variable cualitativa que toma k valores, debemos crear k-1 nuevas variables numéricas (de un tipo muy concreto) en su lugar.

Sí tenemos una variable cualitativa que toma 4 categorías (cualitativa nominal)

datos_preproc |> count(region)
# A tibble: 4 × 2
  region        n
  <chr>     <int>
1 Midwest    1021
2 Northeast   210
3 South      1368
4 West        448

Para poder incluirla en la regresión la recodificaremos con la técnica conocida como one-hot encoding: creamos nuevas variables binarias (si tenemos \(k=4\) modalidades, crearemos \(k-1 = 3\) variables) que nos diga si el registro toma dicho valor o no (hacemos una menos para evitar colinealidad)

# Escribe el código que consideres
datos_preproc2 <-
  datos_preproc |> 
  mutate(region_Northeast = region == "Northeast",
         region_south = region == "South",
         region_West = region == "West") |> 
  select(-region)

Esto lo podemos hacer de manera automática con dummy_cols del paquete {fastDummies}

# Escribe el código que consideres
datos_preproc2 <-
  datos_preproc |>
  fastDummies::dummy_cols(select_columns = "region",
                          remove_first_dummy = TRUE) |> 
  select(-region)

Ejercicio 3 (1.5 puntos)

Ejecuta el código que consideres para realizar una selección inicial de variables de manera que se mitiguen los efectos adversos de la colinealidad. Comenta todo lo que consideres sobre. Tras ello separa las observaciones en dos datasets distintos: un dataset train y un segundo dataset test. Para ello realiza un muestreo aleatorio de manera que train contenga el 80% de los datos y test el 20% restante. IMPORTANTE: no cambies la semilla.

Vamos a realizar un análisis de correlaciones previo

datos_preproc2 |>
  corrr::correlate()
Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'
# A tibble: 15 × 16
   term      deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
   <chr>         <dbl>     <dbl>      <dbl>          <dbl>       <dbl>     <dbl>
 1 deathRate  NA         -0.429     -0.120          0.429     -0.0223   -0.00426
 2 medIncome  -0.429     NA          0.236         -0.789      0.0440   -0.117  
 3 popEst20…  -0.120      0.236     NA             -0.0653     0.0557   -0.176  
 4 povertyP…   0.429     -0.789     -0.0653        NA         -0.0557   -0.193  
 5 studyPer…  -0.0223     0.0440     0.0557        -0.0557    NA        -0.0317 
 6 MedianAge  -0.00426   -0.117     -0.176         -0.193     -0.0317   NA      
 7 AvgHouse…  -0.0369     0.112      0.110          0.0743    -0.00407  -0.356  
 8 PercentM…  -0.267      0.355     -0.160         -0.643     -0.0381    0.429  
 9 PctHS25_…   0.405     -0.471     -0.312          0.194     -0.0851    0.330  
10 PctUnemp…   0.378     -0.453      0.0508         0.655     -0.0320   -0.127  
11 PctPriva…  -0.386      0.724      0.0527        -0.823      0.0925    0.0687 
12 BirthRate  -0.0874    -0.0102    -0.0577        -0.0123     0.0107   -0.0999 
13 region_N…  -0.0646     0.195      0.135         -0.152      0.0143    0.0621 
14 region_S…   0.331     -0.300     -0.0548         0.421     -0.108    -0.118  
15 region_W…  -0.279      0.114      0.0897        -0.0530    -0.00804  -0.0405 
# ℹ 9 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
#   PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
#   BirthRate <dbl>, region_Northeast <dbl>, region_South <dbl>,
#   region_West <dbl>
datos_preproc2 |>
  cor() |> 
  corrplot::corrplot(method = "color")

Concretamente vamos a fijarnos en las correlaciones vs la variable objetivo

datos_preproc2 |>
  corrr::correlate() |> 
  corrr::focus("deathRate") |> 
  filter(abs(deathRate) < 0.05)
Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'
# A tibble: 3 × 2
  term             deathRate
  <chr>                <dbl>
1 studyPerCap       -0.0223 
2 MedianAge         -0.00426
3 AvgHouseholdSize  -0.0369 

Como vemos tenemos 3 variables con una correlación con la variable objetivo por debajo de 0.05 en valor absoluto (prácticamente incorreladas): studyPerCap, MedianAge y AvgHouseholdSize. Podríamos no haberlo hecho (ya que no se pedía el enunciado pero si lo has hecho, es lo correcto) o incluso ser más exigente y subir el umbral de correlación. Vamos a quitar de momento solo esas 3 que presentan una aparente incorrelación lineal (!= independencia) con la objetivo

datos_cor <-
  datos_preproc2 |>
  select(-studyPerCap, -MedianAge, -AvgHouseholdSize)
datos_cor |> 
  cor() |> 
  corrplot::corrplot(method = "color")

Tras ello vemos que hay predictoras altamente correlacionadas entre sí:

  • medIncome con povertyPercent con $-0.789 $ (lógico: más ingresos, menos pobre)

  • PctPrivateCoverage y medIncome con \(0.724\)

  • PctPrivateCoverage y povertyPercent con \(-0.822\)

Para decidir cuál quitamos vamos a basarnos en dos criterios: aquella cuya relación lineal con el conjunto de las demás sea más alta y que tenga menos dependencia con la variable objetivo (tienen más o menos similar entre \(0.43\) y \(0.38\) en valor absoluto)

Este análisis de correlaciones solo nos permite ver relaciones 2 a 2 entre las variables, pero no dependencias lineales más complejas, y para eso vamos a realizar un análisis de multicolinealidad con check_collinearity()

ajuste_saturado <- lm(data = datos_cor, formula = deathRate ~ .)
performance::check_collinearity(ajuste_saturado)
# Check for Multicollinearity

Low Correlation

                 Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
            medIncome 4.22 [3.97, 4.50]         2.06      0.24     [0.22, 0.25]
           popEst2015 1.20 [1.15, 1.25]         1.09      0.84     [0.80, 0.87]
       PercentMarried 2.30 [2.18, 2.44]         1.52      0.43     [0.41, 0.46]
         PctHS25_Over 1.72 [1.64, 1.81]         1.31      0.58     [0.55, 0.61]
 PctUnemployed16_Over 2.06 [1.95, 2.17]         1.43      0.49     [0.46, 0.51]
   PctPrivateCoverage 4.20 [3.94, 4.47]         2.05      0.24     [0.22, 0.25]
            BirthRate 1.07 [1.04, 1.13]         1.04      0.93     [0.89, 0.96]
     region_Northeast 1.27 [1.22, 1.33]         1.12      0.79     [0.75, 0.82]
         region_South 1.83 [1.74, 1.93]         1.35      0.55     [0.52, 0.58]
          region_West 1.60 [1.53, 1.68]         1.27      0.62     [0.59, 0.65]

Moderate Correlation

           Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
 povertyPercent 6.83 [6.40, 7.29]         2.61      0.15     [0.14, 0.16]

Parece que de lsa 3 la que “colisiona” más con el resto de predictoras es povertyPercent por lo que procederemos retirarla

datos_cor <-
  datos_cor |> 
  select(-povertyPercent)
ajuste_saturado <- lm(data = datos_cor, formula = deathRate ~ .)
performance::check_collinearity(ajuste_saturado)
# Check for Multicollinearity

Low Correlation

                 Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
            medIncome 3.05 [2.88, 3.24]         1.75      0.33     [0.31, 0.35]
           popEst2015 1.20 [1.15, 1.25]         1.09      0.84     [0.80, 0.87]
       PercentMarried 1.76 [1.67, 1.85]         1.33      0.57     [0.54, 0.60]
         PctHS25_Over 1.69 [1.61, 1.78]         1.30      0.59     [0.56, 0.62]
 PctUnemployed16_Over 1.98 [1.89, 2.10]         1.41      0.50     [0.48, 0.53]
   PctPrivateCoverage 3.48 [3.28, 3.70]         1.87      0.29     [0.27, 0.31]
            BirthRate 1.07 [1.04, 1.13]         1.04      0.93     [0.89, 0.96]
     region_Northeast 1.26 [1.21, 1.32]         1.12      0.79     [0.76, 0.83]
         region_South 1.82 [1.73, 1.91]         1.35      0.55     [0.52, 0.58]
          region_West 1.59 [1.52, 1.67]         1.26      0.63     [0.60, 0.66]

Por último vamos a realizar la partición pedida

# Escribe el código que consideres
set.seed(12345)
datos_split <- initial_split(datos_cor, prop = 0.8)
train <- training(datos_split)
test <- testing(datos_split)

Ejercicio 4 (1 punto)

No usaremos el dataset test hasta el final para evaluar las predicciones (con un dataset que el modelo no conoce). Con el dataset train ejecuta el ajuste de regresión lineal multivariante saturado (habiendo hecho lo pedido en los ejercicios anteriores) y comenta de manera MUY RESUMIDA la salida (SOLO lo que puedas interpretar hasta ahora, y al grano, que se te acaba el tiempo)

# Escribe el código que consideres
ajuste_saturado <- lm(data = train, formula = deathRate ~ .)
ajuste_saturado |> summary()

Call:
lm(formula = deathRate ~ ., data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-103.288  -11.893    0.604   11.863  158.564 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.605e+02  7.315e+00  21.941  < 2e-16 ***
medIncome            -1.690e-04  6.480e-05  -2.608  0.00917 ** 
popEst2015           -1.060e-06  1.352e-06  -0.784  0.43314    
PercentMarried       -4.777e-01  8.531e-02  -5.599 2.40e-08 ***
PctHS25_Over          1.228e+00  8.342e-02  14.721  < 2e-16 ***
PctUnemployed16_Over  1.669e+00  1.834e-01   9.096  < 2e-16 ***
PctPrivateCoverage   -5.346e-02  7.747e-02  -0.690  0.49016    
BirthRate            -5.073e-01  2.342e-01  -2.167  0.03036 *  
region_Northeast     -3.588e+00  1.956e+00  -1.835  0.06667 .  
region_South          7.577e+00  1.216e+00   6.229 5.51e-10 ***
region_West          -1.168e+01  1.602e+00  -7.290 4.18e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.99 on 2426 degrees of freedom
Multiple R-squared:  0.3736,    Adjusted R-squared:  0.371 
F-statistic: 144.7 on 10 and 2426 DF,  p-value: < 2.2e-16

Algunos comentarios que se podrían hacer:

  • El \(R^2\) es 0.3736 y el ajustado \(0.371\), lo que nos da un indicio de que tampoco tenemos un sobreajuste muy grande (tras quitar lo que hemos quitado ya, claro)

  • Ese \(R^2\) prácticamente duplica el univariante de la primera entrega

  • A mayores ingresos, disminuye la tasa de mortalidad: por cada 10000$ que suba los ingresos medianos, la tasa de mortalidad baja 1.69 puntos.

  • La educación y empleo importa tanto que por cada % que se incrementa en PctHS25_Over (es decir, más gente cuyo máximo título es solo bachillerato) o por cada % que se incrementa PctUnemployed16_Over, la tasa de mortalidad sube 1.2 y 1.66 puntos respectivamente (menos formación o menos empleo –> peores condiciones de vida –> más mortalidad)

  • El código postal importa: simplemente por el hecho de estar en el Sur, tu tasa de mortalidad sube ¡7 puntos!

NO SE PUEDE COMENTAR LOS P-VALORES YA QUE NO SABEMOS SI LA DIAGNOSIS ES CORRECTA

Ejercicio 5 (1.75 puntos)

Realiza una selección de modelos en base a los criterios BIC y AIC. Recuerda que para evitar incompatibilidad con otros paquetes, no debes hAcer library(MASS) sino MASS::… Comenta e interpreta todo lo que sepas de las salidas generadas. ¿Cuál penaliza más? ¿Por qué? ¿Qué ventajas tiene la selección de BIC? Interpreta los coeficientes de ambos modelos y comenta las diferencias (en caso de que existiesen)

# Escribe el código que consideres
ajuste_AIC <- MASS::stepAIC(ajuste_saturado, direction = "both",
                            k = 2)
Start:  AIC=15073.65
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + PctPrivateCoverage + BirthRate + region_Northeast + 
    region_South + region_West

                       Df Sum of Sq     RSS   AIC
- PctPrivateCoverage    1       230 1172938 15072
- popEst2015            1       297 1173005 15072
<none>                              1172708 15074
- region_Northeast      1      1627 1174335 15075
- BirthRate             1      2269 1174977 15076
- medIncome             1      3287 1175995 15078
- PercentMarried        1     15155 1187863 15103
- region_South          1     18757 1191465 15110
- region_West           1     25688 1198396 15124
- PctUnemployed16_Over  1     39994 1212702 15153
- PctHS25_Over          1    104760 1277467 15280

Step:  AIC=15072.13
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + BirthRate + region_Northeast + region_South + 
    region_West

                       Df Sum of Sq     RSS   AIC
- popEst2015            1       251 1173189 15071
<none>                              1172938 15072
- region_Northeast      1      1537 1174475 15073
+ PctPrivateCoverage    1       230 1172708 15074
- BirthRate             1      2125 1175064 15074
- medIncome             1      6658 1179596 15084
- PercentMarried        1     15328 1188266 15102
- region_South          1     23971 1196909 15119
- region_West           1     26930 1199868 15125
- PctUnemployed16_Over  1     48377 1221315 15169
- PctHS25_Over          1    104776 1277714 15279

Step:  AIC=15070.65
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + 
    BirthRate + region_Northeast + region_South + region_West

                       Df Sum of Sq     RSS   AIC
<none>                              1173189 15071
- region_Northeast      1      1623 1174812 15072
+ popEst2015            1       251 1172938 15072
+ PctPrivateCoverage    1       184 1173005 15072
- BirthRate             1      2108 1175297 15073
- medIncome             1      7272 1180460 15084
- PercentMarried        1     15080 1188269 15100
- region_South          1     24140 1197329 15118
- region_West           1     26985 1200174 15124
- PctUnemployed16_Over  1     48126 1221315 15167
- PctHS25_Over          1    109109 1282298 15285
ajuste_AIC <-
  lm(data = train,
     formula = 
       deathRate ~ medIncome + PercentMarried + PctHS25_Over +
       PctUnemployed16_Over + BirthRate + region_Northeast +
       region_South + region_West)
ajuste_AIC |> summary()

Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + BirthRate + region_Northeast + region_South + 
    region_West, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-103.450  -11.926    0.524   11.870  158.525 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.572e+02  6.065e+00  25.912  < 2e-16 ***
medIncome            -2.012e-04  5.187e-05  -3.879 0.000108 ***
PercentMarried       -4.713e-01  8.436e-02  -5.587 2.58e-08 ***
PctHS25_Over          1.238e+00  8.236e-02  15.027  < 2e-16 ***
PctUnemployed16_Over  1.705e+00  1.709e-01   9.980  < 2e-16 ***
BirthRate            -4.851e-01  2.323e-01  -2.089 0.036829 *  
region_Northeast     -3.564e+00  1.945e+00  -1.833 0.066953 .  
region_South          7.926e+00  1.121e+00   7.068 2.04e-12 ***
region_West          -1.133e+01  1.517e+00  -7.473 1.09e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.98 on 2428 degrees of freedom
Multiple R-squared:  0.3733,    Adjusted R-squared:  0.3712 
F-statistic: 180.8 on 8 and 2428 DF,  p-value: < 2.2e-16
ajuste_BIC <- MASS::stepAIC(ajuste_saturado, direction = "both",
                            k = log(nrow(train)))
Start:  AIC=15137.43
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + PctPrivateCoverage + BirthRate + region_Northeast + 
    region_South + region_West

                       Df Sum of Sq     RSS   AIC
- PctPrivateCoverage    1       230 1172938 15130
- popEst2015            1       297 1173005 15130
- region_Northeast      1      1627 1174335 15133
- BirthRate             1      2269 1174977 15134
- medIncome             1      3287 1175995 15136
<none>                              1172708 15137
- PercentMarried        1     15155 1187863 15161
- region_South          1     18757 1191465 15168
- region_West           1     25688 1198396 15182
- PctUnemployed16_Over  1     39994 1212702 15211
- PctHS25_Over          1    104760 1277467 15338

Step:  AIC=15130.11
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + BirthRate + region_Northeast + region_South + 
    region_West

                       Df Sum of Sq     RSS   AIC
- popEst2015            1       251 1173189 15123
- region_Northeast      1      1537 1174475 15126
- BirthRate             1      2125 1175064 15127
<none>                              1172938 15130
- medIncome             1      6658 1179596 15136
+ PctPrivateCoverage    1       230 1172708 15137
- PercentMarried        1     15328 1188266 15154
- region_South          1     23971 1196909 15172
- region_West           1     26930 1199868 15178
- PctUnemployed16_Over  1     48377 1221315 15221
- PctHS25_Over          1    104776 1277714 15331

Step:  AIC=15122.84
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + 
    BirthRate + region_Northeast + region_South + region_West

                       Df Sum of Sq     RSS   AIC
- region_Northeast      1      1623 1174812 15118
- BirthRate             1      2108 1175297 15119
<none>                              1173189 15123
- medIncome             1      7272 1180460 15130
+ popEst2015            1       251 1172938 15130
+ PctPrivateCoverage    1       184 1173005 15130
- PercentMarried        1     15080 1188269 15146
- region_South          1     24140 1197329 15165
- region_West           1     26985 1200174 15170
- PctUnemployed16_Over  1     48126 1221315 15213
- PctHS25_Over          1    109109 1282298 15332

Step:  AIC=15118.41
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + 
    BirthRate + region_South + region_West

                       Df Sum of Sq     RSS   AIC
- BirthRate             1      1742 1176554 15114
<none>                              1174812 15118
+ region_Northeast      1      1623 1173189 15123
+ popEst2015            1       337 1174475 15126
+ PctPrivateCoverage    1        98 1174714 15126
- medIncome             1      9521 1184333 15130
- PercentMarried        1     13708 1188520 15139
- region_West           1     25366 1200178 15163
- region_South          1     31422 1206234 15175
- PctUnemployed16_Over  1     47178 1221990 15207
- PctHS25_Over          1    107638 1282450 15324

Step:  AIC=15114.22
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + 
    region_South + region_West

                       Df Sum of Sq     RSS   AIC
<none>                              1176554 15114
+ BirthRate             1      1742 1174812 15118
+ region_Northeast      1      1257 1175297 15119
+ popEst2015            1       308 1176246 15121
+ PctPrivateCoverage    1        26 1176528 15122
- medIncome             1      8811 1185365 15125
- PercentMarried        1     15444 1191998 15138
- region_West           1     25515 1202069 15159
- region_South          1     32512 1209067 15173
- PctUnemployed16_Over  1     47741 1224296 15203
- PctHS25_Over          1    109486 1286040 15323
ajuste_BIC <-
  lm(data = train,
     formula = 
       deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + region_South + region_West)
ajuste_BIC |> summary()

Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + region_South + region_West, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-105.425  -12.091    0.501   11.766  159.447 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.541e+02  5.944e+00  25.933  < 2e-16 ***
medIncome            -2.141e-04  5.020e-05  -4.266 2.07e-05 ***
PercentMarried       -4.624e-01  8.187e-02  -5.648 1.82e-08 ***
PctHS25_Over          1.233e+00  8.202e-02  15.038  < 2e-16 ***
PctUnemployed16_Over  1.694e+00  1.706e-01   9.930  < 2e-16 ***
region_South          8.706e+00  1.062e+00   8.194 4.02e-16 ***
region_West          -1.072e+01  1.476e+00  -7.259 5.21e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 22 on 2430 degrees of freedom
Multiple R-squared:  0.3715,    Adjusted R-squared:   0.37 
F-statistic: 239.4 on 6 and 2430 DF,  p-value: < 2.2e-16
performance::compare_performance(ajuste_saturado, ajuste_AIC, ajuste_BIC)
# Comparison of Model Performance Indices

Name            | Model |   AIC (weights) |  AICc (weights) |   BIC (weights) |    R2 | R2 (adj.) |   RMSE |  Sigma
-------------------------------------------------------------------------------------------------------------------
ajuste_saturado |    lm | 21991.6 (0.154) | 21991.7 (0.151) | 22061.1 (<.001) | 0.374 |     0.371 | 21.936 | 21.986
ajuste_AIC      |    lm | 21988.6 (0.690) | 21988.6 (0.691) | 22046.5 (0.013) | 0.373 |     0.371 | 21.941 | 21.982
ajuste_BIC      |    lm | 21991.5 (0.156) | 21991.6 (0.158) | 22037.9 (0.987) | 0.372 |     0.370 | 21.972 | 22.004
  • Ajuste AIC: ha eliminado las variables popEst2015 y PctPrivateCoverage, empeorando ligeramente \(R^2\) pero incrementando ligeramente el \(R^2\) (menos variables pero más necesarias)

  • Ajuste BIC: ha eliminado además las variables BirthRate y region_Northeast, empeorando ligeramente \(R^2\)

La diferencia de ambos es que el segundo es más “agresivo” penalizando. ¿La ventaja? Que asintóticamente (si \(n\) fuese incrementándose) tenemos garantizado que el BIC acaba eligiendo el modelo real que subyace en los datos (consistencia)

Al final tenemos 3 modelos con prácticamente la misma cantidad de información explicada, pero uno tiene 10 variables, otro 8 y otro 6 variables (con casi la misma calidad, pero con casi la mitad de predictoras) –> nos qedamos con el ajuste BIC para la diagnosis.

Ejercicio 6 (1.5 puntos)

Compara la calidad de los 3 modelos (saturado, BIC y AIC) en train. Quédate con uno (justifica) y comprueba si se cumplen las hipotesis necesarias (fase de diagnosis SOLO en ese elegido). ¿Se cumplen las hipótesis? Argumenta porqué, más allá de lo que valga luego la bondad de ajuste, es importante que se cumplan.

# Escribe el código que consideres
performance::check_model(ajuste_BIC)

En un primer vistazo visual vemos que es bastante probable que muchas hipótesis se cumplan o estén derca de complirse (teniendo en cuenta que no hemos tratado los outliers que pueden estar influenciando)

Linealidad

ggplot(tibble("fitted" = ajuste_BIC$fitted.values,
              "residuals" = ajuste_BIC$residuals),
       aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

A la izquierda y arriba se observan algunas observaciones que probablemente sean outliers peor el resto si parece presentar un patrón independiente entre residuos y valores ajustados, no viéndose ninguna tendencia o modelo posible a ajustar entre ellos, así que podríamos dar por válida la linealidad

Homocedasticidad

performance::check_heteroscedasticity(ajuste_BIC)
Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).
ggplot(tibble("id" = 1:length(ajuste_BIC$residuals),
              "residuals" = ajuste_BIC$residuals),
       aes(x = id, y = residuals)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

A pesar de que el p-valor nos dice que se rechaz ala homocedasticidad, los errores no presentan ningún patrón ni aumentan su varianza, estando prácticamente todos dentro de una banda constante (salvo unos que probablemente sean outliers y deberíamos analizarlos de manera particular)

Normalidad

performance::check_normality(ajuste_BIC)
Warning: Non-normality of residuals detected (p < .001).
ggplot(tibble("residuals" = ajuste_BIC$residuals)) +
  geom_qq(aes(sample = residuals)) +
  geom_qq_line(aes(sample = residuals)) +
  theme_minimal()

Normalidad no parece cumplirse. Recuerda que aquí influye:

  • tamaño muestral
  • outliers que hace que en los extremos de los quantiles la distribució se aleje
ggplot(tibble("residuals" = ajuste_BIC$residuals,
              "normal" = rnorm(n = length(ajuste_BIC$residuals), mean = 0,
                               sd = sd(ajuste_BIC$residuals)))) +
  geom_density(aes(x = residuals), color = "blue") +
  geom_density(aes(x = normal), color = "red") +
  theme_minimal()

Vemos como la distribución de los residuos es algo más apuntada que lo que correspondería a una distribución normal (probablemente por los outliers que contienen nuestros datos)

Correlación

performance::check_autocorrelation(ajuste_BIC)
OK: Residuals appear to be independent and not autocorrelated (p = 0.612).
acf(ajuste_BIC$residuals)

Vemos tanto numérica como visualmente uqe los residuos están incorrelados (con acf() debemos obtener lo contrario a lo que tendríamos al inicio de un análisis de series temporales: los residuos retardados sin correlacion entre ellos)

Predicciones simulación

performance::check_predictions(ajuste_BIC)
Warning: Maximum value of original data is not included in the
  replicated data.
  Model may not capture the variation of the data.

Las simulaciones (bajo la hipotesis de que el modelo ajustado fuese el modelo real subyacente) están ligeramente por debajo del dato observado, probablemente producido por ese apuntamiento que se observaba en los residuales estimados.

Ejercicio 7 (0.75 puntos)

Asume que se cumplen las hipótesis (en caso de que no se cumpliesen) e interpreta todo lo que sepas sobre la segunda tabla de la salida de la regresión (la de inferencia de los parámetros). Elimina variables si su efecto no fuese significativo.

ajuste_BIC |> summary()

Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + region_South + region_West, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-105.425  -12.091    0.501   11.766  159.447 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.541e+02  5.944e+00  25.933  < 2e-16 ***
medIncome            -2.141e-04  5.020e-05  -4.266 2.07e-05 ***
PercentMarried       -4.624e-01  8.187e-02  -5.648 1.82e-08 ***
PctHS25_Over          1.233e+00  8.202e-02  15.038  < 2e-16 ***
PctUnemployed16_Over  1.694e+00  1.706e-01   9.930  < 2e-16 ***
region_South          8.706e+00  1.062e+00   8.194 4.02e-16 ***
region_West          -1.072e+01  1.476e+00  -7.259 5.21e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 22 on 2430 degrees of freedom
Multiple R-squared:  0.3715,    Adjusted R-squared:   0.37 
F-statistic: 239.4 on 6 and 2430 DF,  p-value: < 2.2e-16

Asumiendo que las hipótesis son ciertas, interpretando los p-valores, vemos todas las covariables presentan un efecto lineal significativo así como a nivel global según el ANOVA realizado

ajuste_BIC |> anova()
Analysis of Variance Table

Response: deathRate
                       Df  Sum Sq Mean Sq F value    Pr(>F)    
medIncome               1  362519  362519 748.729 < 2.2e-16 ***
PercentMarried          1   26011   26011  53.722 3.128e-13 ***
PctHS25_Over            1  152438  152438 314.839 < 2.2e-16 ***
PctUnemployed16_Over    1   53142   53142 109.757 < 2.2e-16 ***
region_South            1   75831   75831 156.618 < 2.2e-16 ***
region_West             1   25515   25515  52.697 5.214e-13 ***
Residuals            2430 1176554     484                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Así lo confirmamos con el anova secuencial: en cada paso, se calcula la diferencia de SSR entre el modelo con las variables anteriroes y el modelo con todas más la variable correspondiente

Ejercicio 8 (1.5 puntos)

Por último, utiliza el dataset test para, aplicando los 3 modelos entrenados (saturado, AIC y BIC), predecir las observaciones de test (con cada uno de los modelos. Con las observaciones de test construye un nuevo dataset de 3 columnas: la variable objetivo, la predicción y el modelo usado (saturado, BIC, AIC). Calcula el \(R^2\) en el dataset de test para cada modelo. ¿Cuál funciona mejor en test? ¿Qué % de información explica cada modelo? ¿Cuánto vale su varianza residual?

# Escribe el código que consideres
pred <-
  tibble("modelo" = "saturado",
         "real" = test$deathRate,
         "pred" = predict(ajuste_saturado, test)) |> 
  bind_rows(tibble("modelo" = "AIC",
                   "real" = test$deathRate,
         "pred" = predict(ajuste_AIC, test))) |> 
  bind_rows(tibble("modelo" = "BIC",
                   "real" = test$deathRate,
         "pred" = predict(ajuste_BIC, test)))
pred |> 
  summarise(SSR = sum((real - pred)^2),
            SST = sum((real - mean(real))^2),
            var_res = SSR / nrow(train),
            R2 = 1 - SSR/SST, .by = modelo)
# A tibble: 3 × 5
  modelo       SSR     SST var_res    R2
  <chr>      <dbl>   <dbl>   <dbl> <dbl>
1 saturado 323375. 473810.    133. 0.317
2 AIC      324694. 473810.    133. 0.315
3 BIC      328353. 473810.    135. 0.307

Vemos como el modelo con un \(R^2\) más bajo es el ajuste BIC, pero con muy poca diferencia (algo que ya intuíamos ya que sucedía así en train). En los 3 casos los modelos explican aprox el 30% de información contenida en la variable objetivo